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1. INTRODUCTION 



The determination of correlation functions is an important task both from a theoretical 
and from an experimental point of view. In scattering experiments correlation functions are 
measured, and detailed insight into the microscopic structure of the sample can be obtained. 
Scattering experiments may thus well decide about the adequateness of a given (simplified) 
theoretical model, once its correlations functions are calculated. On the other hand, taken for 
granted that a theoretical model adequately describes a certain sample, the knowledge of its 
correlation functions is useful to estimate the quality of the experimental setup. The principal 
problems on the theoretical side are threefold: To deal with many-body interactions between 
the particles, to account for finite temperatures and to carry out the thermodynamic limit of 
a macroscopic sample size. With regard to these intricate difficulties, techniques which allow 
for the exact calculation of correlation functions of at least some typical interacting quantum 
systems in the thermodynamic limit are highly welcome. 

Progress into this direction has recently been achieved by one of the authors in collaboration 
with A. Klumper and A. Seel (GKS) 

BSS. 

In these works, the one-dimensional spin-1/2 
XXZ Heisenberg chain with periodic boundary conditions in a magnetic field, 

H = JJ2 \s*S* +1 + S]^ + A ( S !Sj +1 -\)] -hj^sh (!) 

3=1 L V /J j=l 

was considered. Formulae for the correlation functions were derived by combining the quantum- 
transfer-matrix (QTM) approach (J, 0, If], [?1 to the thermodynamics of integrable systems with 
certain results for matrix elements of the XXZ-ch&m N, M which proved to be useful before 

n n 

jlOL II 1| at T = 0. In the QTM approach, the free energy of the system is expressed in 
terms of an auxiliary function defined in the complex plane []]]. This auxiliary function is 
determined uniquely as a solution of a non-linear integral equation. GKS found that essentially 
the same auxiliary function can be used to calculate m-point correlation functions for arbitrary 
temperatures T, again in the thermodynamic limit. As a result, all the correlation functions 
of spatial range m are given in terms of m-fold integrals over combinations of the auxiliary 
function taken at different integration variables. 

In this work, the integral equations are solved numerically for all temperatures at h — and 
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the multiple integrals are calculated for m — 2, 3, where we restrict ourselves to A > 1 for m — 3. 
In particular, we focus on (5'JS'J +m _ 1 )T and on the emptiness formation probability P(m). In 
order to estimate the quality of our data we compare them with other exact results: first 
of all the nearest-neighbor correlators can be obtained independently of the multiple integral 
representation by taking the derivative of the free energy with respect to A. Furthermore, 
very recently [l2| a high temperature expansion of the emptiness formation probability -P(3) 
at A = 1 has been performed up to the order (J/T) 42 , starting from the multiple integral 

nnhnnn 

representation. At T = 0, closed expressions are available |13j, |1J, UM UM, 111 UM as well. There 
are, moreover, the XX-limit (A — > 0) and the Ising limit ( J — > 0, A — > ock J A fixed) where 

nnn 

the correlation functions are known over the whole range of temperatures [19|, |20|, |21[ . All these 
independent results provide useful tests for the accuracy of our numerics. The main error we 
observe is due to the discretization of the integrals, an error which, in principle, can be made 
arbitrarily small by increasing the numerical accuracy. 

Besides the expected crossover behaviour from low to high temperatures, we find a surprising 
feature in (>5J5'J +m _ 1 )T for A > 1: the antiferromagnetic correlation is maximal at an interme- 
diate temperature To, and not, as naively expected, at T = 0. We explain this behaviour from 
a competition between quantum and thermodynamical fluctuations. 

This paper is organized as follows: in the next section we give some details about the 
numerical procedure. Results are presented separately in the third section. The article ends 
with a summary and an outlook. 

2. NUMERICAL PROCEDURE 

Let Ai... m be an operator which acts on sites 1, . . . , m of the spin chain. In order to calculate 
{A-i...m)Ti the operator A\_ m is expanded in terms of the elementary matrices e^, 






e 1 



5 °2 — ) °2 




such that 

A _ APx-Pm p ai a 2 a, 

Sil...m — /1 a 1 ...a m e l / g 1 e 2/3 2 • • ■ e m/3. 
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Then it suffices to calculate the expectation value (e 1 ^ . . .e m ^™)r which defines the density 
matrix of the chain segment consisting of sites 1, . . . , m. 



In |3| the following multiple integral representation for the density matrix was has been 
obtained 3^|: 



)a+\ 

n 

j=i 



do;, 



27ri(l + a(iVj)) 



n 

j=|a+|+l 



do;. 



c 2m(l + a(ojj)) 



det r 



x- 



x 



X 



1 d b 



(fe-1)! d£ b 



-G(u a ,£) \ (= 



\a+\ 



JJsinh^ \uj la+hj+l 

3=1 



rj) sinh r 



a+|-i+l 



m— |a~! 



sinh^ "Vla+l+j + V) sinh m " 



\a+\+j 



(2) 



Here (a n )™ =1 and (/3 n )™ =1 denote sequences of up- and down-spins, where 1 =f and 2 =|. The 
position of the jth up-spin (down-spin) in the sequence (a n )™ =1 ((/? n )n=i) i s a t W7) an d the 
number of up- (down-) spins in the corresponding sequences is defined as |a + | The 
contour C in © is given in it depends on the value of A =: cosh 77. The parameter rj may 
be purely real and positive (massive case, A > 1), or imaginary with i] =: i^y, ir/2 > 7 > 
(massless case, < A < 1). The auxiliary functions a, a, G which occur in the integrand in (J2J) 
are solutions of the following integral equations: 



In o(A) = -(3h - —— kv (A + 77/2) 



2vri 



7 K r? (A — to) ln(l + a(o>)) 



G(A,0 = %(A-e-^/2)+ / ^k v (\-uj)^p 
2 J c 2iri 1 + a 

sinh(2?7) 



G(cu,Q_ 



(3) 
(4) 



a : = 



sinh(A + 77) sinh(A — 77) 
l/o, 



with /3 := 1/T. Note that the T- and /Vdependence of the correlation functions enters through 
the auxiliary functions, whereas the m-dependence is encoded in the m-fold integral. Thus the 
evaluation of (j2J) requires two steps: solving equations 0, (@J) and then calculating the multiple 
integrals. In this work, the following correlation functions will be calculated numerically at 
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FIG. 1: Integration contours a) in the massive b) in the massless case. 



h = 0: 

P(2) := (e il 1 e i+ ii>r = ^+(^ + i>T 

P(3) := (e^ Je i+1 ie i+2 i) T = - + -(SjS? +1 ) T + (SjSj +2 ) T 
(S*S* +2 ) T = \{e 3 l{e ]+l \ + e J+l l)e 0+2 \) T . 

The function P(m) = (e 1 \ . . . e m \}r is referred to as the emptiness formation probability. 

For the numerical treatment, it is convenient to reformulate equations (j3J), © such that the 
integrations are done along strai ght lines parallel to the real or imaginary axes. This procedure 



is due to Klumper (cf. the book |22j and references therein). We review the main steps here in 
order to introduce a notation which will prove to be convenient for our purposes. 

The integration contour C is chosen according to fig. ^Jl) ^p)) in the massive (massless) 
case. In the massive (massless) regime, Re C± = ±r//2 =F e (Im C± = ±7/2 =1= e), where e is a 
small positive quantity; the limit e — > is included implicitly. According to eqs. (@J) the 
functions In a, G are given in the part of the complex plane enclosed by C, once they have been 
calculated on C. However, for the calculation of correlators, (J2J, these functions are only needed 
on C. Since in the massive case, all functions are periodic with period m, integrals along Ci ; 2 in 
opposite directions cancel each other. We thus have to calculate the auxiliary functions along 
two straight lines C± parallel to the imaginary (real) axis in the massive (massless) case. 



Let us first focus on the calculation of the auxiliary functions in the massless case, the 
massive case is treated analogously afterwards. Subsequently, we comment on the calculation 
of the multiple integrals with the help of the auxiliary functions. 



2.1. Massless case 

We first concentrate on0<7<7r/2 and deal with the isotropic case 7 = at the end of 
this section. It is convenient to define 

b(x) := a(x + vy/2), b(x) := a{x-vy/2) (5a) 
03 := 1 + b, 03:= 1 + b . (5b) 

Note that b = b*\h^-h- We now perform a "particle- hole-transformation" in Q by substituting 

ln(l + a(x + i7/2)) = ln03(x) on C + 

— (6) 
ln(l + a(x - vy/2)) = In 03 (a:) + lnb(x - i7) on C_ . 

Taking account of the simple pole of In b(x) at x — —vy/2, we arrive at 



In b(x) = -/3h + J P sm ^^ k h/2 {x) + k h * In 03 - k<J ] * In 23 - k h * In b (x) (7) 



with the definitions 



f^\x) := k h (x + ie) 



oc 



[f*g](x) ■= / f(x-y)g(y)dy. (8) 

The integrations in (J7J) are now convolutions along the real axis, such that one can express In b 
in terms of In 03, In 03 by solving an algebraic equation in Fourier space and transforming back 
to direct space. The result is 

In b(x) = , - . (3h - J/3sm7 d(x) + \k * In 03 - «+ * In 03 1 (x) (9a) 

2(ir — 7) 2 L J 

lnb = lnb*| h ^ (9b) 



(3 



with 



d(x) 
k(x) 

K+(x) 



n 



7 cosh -x 

sinh (| -7) ke ikx dk 
, 2cosh|fc sinh (f - fc 2^ 



(10) 



k(x + i7 — ie) 



Note that the shift by e in (J5ajl is necessary to ensure integrability. In a similar fashion, eq. (@J) 
is treated. Therefore we introduce the functions 



G b (x) := G{x + i 7 /2, , GV := -G(x - i 7 /2, 0, 



'li: 



where we suppressed the dependence on £ on the left-hand sides. It is now convenient to 
substitute 

G-r 



G 



1 + a 



1 + b 

Gh 

1 + b- 1 

Gir 



on C_ 



Gu on C A 



(12) 



G 



l + o 



1 + b 

Gu 



1 + b" 



3- — Gf on C_ 



on C+ . 



(13) 



Then 



G[,(x) = id(x — £) 
1 



K * 



Gh 



1 + b' 



G, 



7 — 1 
1 + b 



[x 



(14a) 
(14b) 



Using the substitution rules ()12|1. (|13jh the general density matrix element (J2J is expressed 
in terms of G(,/(l + b" 1 ), G^/(l + b ), G(,, G^. This expression is rather lengthy: For an 
m-point-function, one obtains a sum of 3 m -many terms, each one an m-fold integral. Here we 
would like to spare the reader the general expression, but rather make the following comments: 



Eqs. (JyJ), ()14a|) are solved numerically by iteration; the convolutions are done by apply- 
ing the Fast- Fourier- Transform algorithm. Errors then arise from the truncation and the 
discretization of the integrals. 



• In (J2J) derivatives of G bb with respect to £ enter. These can be calculated straightforwardly by 
taking the n-th derivative of ()14|) with respect to £ which results in linear integral equations 
for d-G b - b . 

• The formulation in terms of the functions G bb is particularly suited for low temperatures. At 
T = h = 0, 1/(1 + b -1 ) = 1/(1 + b 1 ) = such that G b = G b is given by the inhomogeneity 
in ()14a|) . Furthermore, only one of the 3 m -many terms in © remains (namely the one 
containing only combinations of Gh, G-r). Then in this term all functions are known. This 

n n 

is the well-known result of Jimbo et al. 1241 which has first been obtained for A > 1 in 23 

n 

(for a proof in the critical case and an extension to h ^ see also [10]). 

The functions for the isotropic case 7 = are obtained by rescaling x — > jx, G bb — > r yG bb . 
Then in (J2J) algebraic functions occur instead of hyperbolic ones, and the integration kernel k 
takes the form 



K(X) 



dk 



1 + el fc l 27r' 



2.2. Massive case 

Let us go through the changes which have to be done in the massive case with respect to 
the massless case. Analogously to (15 all. (|5bl). we define the functions 



b(x) = a(ix + rj/2), b(x) = a(ix — rj/2) 
<B := 1 + b, »:= 1 + b . 

We use the same symbols b, b etc. as in the massless case. Note that here —tt/2 < x < it/2 and 
the b-functions are periodic on the real axis with period it. Accounting for the pole of lnb(x) 
at x = vq/2 we find the following equations 

lnb(x) = -ph - J/3smhr ? k ( x ) + — [ft * ln<B - ft„ * lnb - * ln»1 (x) (15) 

2 2tt 1 

\nb(x) = \nb*\ h ^_ h , 



S 



where the kernel and the convolutions are now defined as 

sinh(2?7) 



sin(x + irj) sin(x — irj) 

/■tt/2 

[f*g](x) = / f(x-y)g(y)dy. (16) 

J —tt/2 

Because of periodicity, eq. (fT5)l is manipulated further in Fourier space after performing a 
discrete Fourier transformation. Upon transforming back one obtains 

lnb(x) = -^-^^Ld(x)+ [«;*ln<B-K_*ln^] (x) (17) 
i(x) = ^dn(^,i^ (18) 



7T V 7T 7T 

oo 



gi2nx 



7r(l + e 2, ?H^ 

n=— oo ' 



K-(x) := k(x — IT) + ie). (19) 

The driving term ()18|) is written in terms of the Jacobi elliptic function dn(x, r) j^j], with 
Fourier series expansion 

n e iTnx/K(r) 

dn(x,r) = — — > ( 2 °) 

2K(t) — ' cos nirr 

v ' n=— oo 



in the strip |3(a;)| < $s(tK(t)) of the complex plane. In (J2(Jj) the constant K(t) is defined 
through K(t) := f $§(0, r) and $3 is one of Jacobi's Thet a- functions j^jj. In very close analogy 
one derives equations for G bb , defined by 

G b (x) := G(ix + V /2,£), G b {x):=-G{ix-r}/2,£), 

where, as in the massless case, the ^-dependence is not noted explicitly in G bb . The symbolical 
substitutions (fT2"j) . (fT3Jl still hold, with C± defined in fig. Hk), so that 



G b (x) = -d{x-£) + 



G b Gj 

K * — K_ * 



(x) (21) 



1 + b- 1 i+r 1 . 

C& — G* b \ h ^-h- (22) 



2.3. Calculating the integrals 

Let us consider eq. where the substitutions implied by (|12jl. (jlHj) have been performed. 
One observes that due to the product of hyperbolic functions in the denominator in eq. ((21), 
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multiple poles may occur. We explain in the following how to treat these poles numerically and 
show that contributions due to these poles cancel each other if there is more than one pole. For 
definiteness we treat the massless case, however, the arguments are directly transferable to the 
massive case. 

The product rii<a<fe<m sinh(cj a - u b - ij) with 

J x a + i 7 /2 on C+ 
uj a = < (23) 
I x a — 17/2 on C- 

acquires a zero if there are j, k such that u>j = Xj + vy/2 — ie, Uk = Xk — ij/2 + ie. This leads 
to a factor 



1 



sinh(a;j — Xk — ie) 



1 



+ m5(xj - x k ) . (24) 



sinh(xj — Xk) 

When taking the principal value numerically, one has to account for the regular part at xj = Xk, 
defined by f reg (xj,Xj) = 0: 

p r ffaxk) d = r>~\ r n^k) dXk _ 2ed f {x }| (25) 

smh( Xj - x k ) J^oo J Xj+e smh(x j - x k ) 

The last terms in eqs. (|24j) . ()25|1 are additional terms caused by the pole. 

Now consider the case of multiple poles. Set Uj = Xj + i 7 /2 in the product 
[\i <a<b<m sinh(u; a — u b — vy) and all other uj k = Xk — i 7 /2. Then if j = m — 1, there is 
one simple pole. If j < m — 1, the product sinh(xj — x m _i — ie) sinh(xj — x m — ie) occurs in the 
denominator. It is not difficult to show that the additional contributions due to these two poles 
either vanish directly or cancel each other. One should be aware of the fact that the 5-function 
in (J2IJ) yields no contribution if / = f reg . This is the case for T = 0, where all functions are 
real and the determinant in vanishes if any two arguments of the auxiliary functions are 
equal. 

We now address the question of possible simplifications of the multiple integral representa- 
tion (J2J). Obviously, in the case m = 2 the double integration has to be equal to one single 
integration, since the nearest-neighbor correlators are obtained by taking the derivative of the 
free energy per lattice site / with respect to the anisotropy A, and the free energy is given by 
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a single integral according to |2(. Especially 



0. 



(SjSj +1 )- J 



1 , d A f 



1 



J 

ii A 

1 u 

(SjSj +1 ) T = (SjSj +1 ) T = Y2 + 3J' 



A 



(26a) 
(26b) 
(26c) 



where u = d^fif) is the inner energy per lattice site and / is given by 

f = e -^[d*\n^ (0). 

Here d(x) is given by eq. (JTUJ) (eq. (fTSj)) in the massless (massive) case; the convolution for both 
cases is defined in eq. (JEJ) (eq. fTHjl ). The ground state energy is given by 

J sin 7 f°° sinh(f-|)A; 



2 

J sinh ?7 



2 cosh 2* sinh ^ 



E 



cosh ?]n 



d/c massless 



massive. 



Eqs. (I26ajl . ()26b|) also allow us to extract the leading order in T at low temperatures. We 
consider here the case of vanishing magnetic field, h = 0. From using dilogarithms, one 
finds for T — > in the massless case 



/ = e -T 2 



7 



3Jsin7 



Thus for T — > 



(SjSj +1 ) T - (S-Sj +1 } T =o + (^7 J) 
{SjSj +1 ) T = (SjSj+i)T=o + (r/ J)' 



2 (1 -7 cot 7) 



3 sin 7 
(7 — cot 7 + 7 cot 2 7) 



6 sin 7 
T 2 

(SjSj +1 ) T = (SjSj +1 ) T = (SjSj +1 ) T =o + g-p, A = 1. 



(27a) 
(27b) 
(27c) 



• 3- 



The correlator ( 1 S , |5 , | +1 )r =0 as a function of A is plotted in fig. 3 of ref. [27|. This result has 
been extended to all temperatures in the range — 1 < A < 1 in fig. 2. Note from (|27aj) - 
(|27cj) that the sign of the T 2 -coefficient is positive, so that the antiferromagnetic correlations 
first decrease with increasing temperature. 
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In the massive case, for finite rj, the leadin 
by a saddle-point integration, very similar to 



; temperature-dependent order may be obtained 
23. Then 



f = -e - A x ' 2 T^ 2 e- B l T 
V 

A = 

2Jk 2 K sinh?] 

T k'K . , 
B = J smh?7 

7T 

with the moduli k 1/2 : = tf 2 (0, irj/n)/ti 3 (0, irj/n) and (k') 1/2 := i? 4 (0, irj/ir)/$ 3 (0, irj/n). It is now 
not difficult to show that 

(S?Sj +1 ) T = (S*S* +1 ) T=0 + CT l /\- B ' T 
{S*S* +1 ) T = (S*S* +1 ) T=0 + DT X I\~ B I T 

C = A'' 2 — ( cothry + V — * -) > 

71 \ cosh 2 7](n - 1/2) J 

D = — (A 1/2 B - JCcoshr]) < . 

At first sight, the last inequality is surprising: The antiferromagnetic correlations of (SjSj +1 )T 
increase with increasing temperature. Since \im.T^oo{Sj Sj +1 )t = 0, we expect a minimum at 
finite temperature. This will be confirmed in the next section. Intuitively, this result can be 
understood as follows: At T = 0, the Ising-like anisotropy leads to a Neel-like ground state, 
which enhances the alignment of the spins along the ^-direction. An initial increase of the 
temperature reduces | (SjSj +1 )T\, so that quantum fluctuations of the spins are favored, and 
| (SjSj +1 )T\ increases. At high enough temperatures the spins become uncorrelated, so that 
both correlation functions decay to zero. 

As far as the low-temperature behavior of correlation functions in the massive regime is 
concerned, we would like to note that using the saddle-point integration technique, one finds 
that 

<e x £ • • • e m yj T = {e x £ . . . e m ->=o + C m T^ 2 e~^ + O (T^e^ , Te~ 2 ^) , 

for arbitrary m, with some coefficient C m . 

Up to now, we were not able to derive eqs. (|26aj) . (j26b|) from (J2J). Still, the latter equations 
provide a useful test for the accuracy of our numerical integration scheme for m = 2, as 
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explained in more detail in the following section. For m > 2, simplifications of (j2J) have not 
been performed yet for any finite temperatures except at A = and A — > oo, see below. 
However, for T = 0, where the integrand is known explicitly, the multiple integrals have been 
simplified considerably in the isotropic regime 17, Q]. In the anisotropic regime, the 
multiple integrals have been reduced to single integrals for m < 4 jisl . Q| . These works allow 
us to estimate the numerical accuracy at T = 0. In the other extreme, at T = oo, the spins are 
uncorrelated so that in this limiting case we also have explicit results which can be checked. 
Note, however, that it is rather non-trivial to reproduce these high-temperature results from 
the general formula (|2*|). This was done only very recently by Tsuboi and Shiroishi in jl^ . 
where a high-temperature expansion of P(m) has been performed to high order in 1/T for the 
isotropic case. This is one further result which will serve as a reference point for estimating the 
accuracy of our numerics for m — 3. 

The cases A = (free fermions) and A — > oo such that J ■ A =: c finite (Ising limit) allow 
for substantial simplifications. Details will be explained in a separate publication 0], here we 
only give the results. For A = 0, that is 7 = tt/2, the convolutions in (|9a|. ()14a|) disappear, so 
that all functions are known explicitly for all temperatures. Then 

"T dp e i(a-%> 



Pirn) = det 



w 2ix 1 + exp (j3J cosp — j3h) 



ah 



in agreement with [20]. In the other extreme, the Ising limit, the driving term and integration 
kernels become constant (d(x) — > 1, k(x) — > l/(27r)) such that the auxiliary functions do not 



depend on x and the integrations are trivial. One finds 

1 + (a*) T 



P(m) 

Al,2 

(Or 



l + (a z ) T 1-(^)t A 2 
2 2 Ai 



m— 1 



(28) 



e-^/Mcosh^i 



sinh 2 ^ — h e^ c 



sinh 



sinh 2 f + e^ c 
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3. RESULTS 



We first present results of the next-neighbor correlation functions, comparing data from the 
double integration with those from the single integration. Results for next-nearest neighbors 
will be given subsequently. Correlation functions for longer distances can be evaluated in the 
same way, however, the numerical cost increases exponentially with the number of integrations. 
That is why for the time being in this feasibility study we do not go beyond m — 3. 

3.1. Nearest neighbors 

Figs. 121 El show P(2), (SjSj +1 )x, as calculated from the double integral, for P(2) together 
with the relative error with respect to the data obtained from the derivative of the free energy 
according to (J26a|) . (j26b|) . The single integration results can be considered as exact; the nu- 
merical error is less than 8 digits for all temperatures. We checked that the relative error for 
{SjSj +l )T vanishes at low temperatures as for P(2); at high temperatures, the absolute error 
is about 10~ 4 . 

In both cases, the error increases with increasing temperature, reaching maximal values at 
T ~ J near the isotropic point. Generally, the error is considerably larger at high than at low 
temperatures. This is due to two reasons: 

• At finite temperature, the auxiliary functions are given numerically from the solutions of 
the integral equations, whereas at T = 0, they are known explicitly. The numerical iteration 
scheme necessarily induces errors due to the truncation and discretization of the integrals. 

• The additional contribution from poles in eq. (|24j) increases with increasing temperature and 
causes an additional error in the integrations. 

We have verified that the error is reduced by decreasing the distance between two successive 
integration points. We conclude that the errors are mainly due to the truncation and discretiza- 
tion of the integrals and can be made arbitrarily small at the expense of increasing computation 
time. 

Considering the low-temperature behavior of {Sj Sj +1 )t, one finds an initial decrease from 
T = with a minimum at T , this is illustrated in fig. El We commented on this phenomenon 
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FIG. 2: a) Emptiness formation probability P(2) for nearest neighbors; the upper curve corresponds 
to the XX limit 77 = i-7r/2. The inset shows the T 2 -coemcient at low temperatures in the massless 
case. Straight lines indicate the exact values according to eqs. (|27a|) . (|27cfl . b) Relative error of P(2). 
The curve labeled by "77 = 5 Ising" shows the relative difference between the rj = 5-values of a) and 
the Ising limit eq. (I2*H|) . 
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In 77/ 

FIG. 3: Correlation function (SjSj +1 ) for nearest neighbors; the lower curve corresponds to r\ = 
i0.87r/2. The inset shows the T 2 -coefficient at low temperatures in the massless case. Straight lines 
indicate the exact values according to eq. ()27b|) . 

in the previous section. From fig. |U] one observes that in the Ising limit, the maximum of 
(S*Sf +1 ) T /(SfSJ +1 ) T=0 is located at T /(AJ) « 0.168. 

3.2. Next-Nearest neighbors 

We calculated next-nearest neighbor correlation functions for the isotropic and the massive 
case. In the massless case, the auxiliary functions in (J2J) are multiplied by a kernel which, de- 
pending on the integration variables, can increase exponentially, so that the auxiliary functions 
have to be determined extremely accurately over the whole integration range. This problem is 
still controllable for nearest neighbors, but becomes severe for m > 2. In the isotropic case, 
however, one deals with algebraic functions, which are much better to handle numerically. At 
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anisotropy parameter A > 1 the integration range is finite, so that those problems do not occur. 

In fig. 0] we show the emptiness formation probability P(3). In the isotropic case, we compare 
our data with the high-temperature expansion of 3] an d with Quantum Monte Carlo (QMC) 
data which were obtained in [2]. The relative error of our data shown in fig. H]is larger than 
in fig. |21 which is due to the fact that the data in (j3J) were calculated by using only half of the 
number of integration points compared to the m = 2 case. We checked at single temperature 
values that the error is of the same order as in the m = 2 case when the calculations are done 
with the same numerical accuracy. However, the computational costs would be unreasonably 
high. At lowest temperatures, a comparison with jl^l I15I Q] yields relative errors of the order 
2 ■ 10- 4 . 

As a second example, {SjSj +2 )T is depicted in fig. At lowest temperatures we again 
compare our data with the explicit results of Q, Q] which yields a relative error of 7 ■ 10~ 3 . 
At highest temperatures, (SjSj +2 )T = 0. The absolute error here is of the order 10 -3 . As in 
the m = 2-case we observe an increase of the antiferromagnetic correlation with increasing T, 
reaching a maximum at some finite temperature T in the massive regime. Fig. |U] illustrates 
this behaviour, where (SjSj_ m+1 )T/(SjSj_ m+1 )T=o is depicted for m = 2, 3 and A > 1. 



4. SUMMARY AND OUTLOOK 

We demonstrated that an accurate numerical evaluation of the formula (j2j) is possible by 
considering the special cases m = 2 with A > and m = 3 with A > 1. The numerical error 
was estimated precisely by comparing with alternative approaches. We have found that this 
error is only due to the discretization and, for rj = ry, the truncation of the integrals. For 
m = 3 in the massless case we encountered problems hard to deal with numerically. Also the 
computational costs of an extension of the numerics to m > 3 become very high if one wants 
to maintain the accuracy. Thus, the next step should be to address the question of in how far 
the multiple integrals can be simplified. Here we would like to mention two approaches. 

In an high-temperature expansion (HTE) has been applied to calculate P{m). The 
excellent agreement of these data with those obtained by numerically evaluating the multiple 
integrals directly (cf. fig. HJ, even at temperatures below the crossover region highlights the 
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FIG. 4: a) Emptiness formation pro bability P(3) for next-nearest neighbors. Shown are results from 



the 3-fold integra 



n 



from QMC [12| and high-temperature expansion (HTE) combined with Pade 
approximation |l2j|. b) Relative error of the data from the 3-fold integral with respect to the QMC 
and the HTE-data. 
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FIG. 5: The correlation function (SJS , | +2 )t f° r V = 0, 1. Note that at low temperatures for A > 1, 
the correlation increases, before it starts to decrease to zero. 
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32| the multiple 



applicability of the HTE method. On the other hand, in the works ^ , 
integrals have been reduced to single integrals at T = 0. The case m = 2 shows that simpli- 
fications are possible also at finite T (compare also Qj), where the integrand in (j2J) is given 
implicitly by integral equations. Further progress may be expected from extending the recent 



results 



for the inhomogeneous generalizations of the correlation functions of the XX Z 



chain to finite T. This issue is the subject of current research. 
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In T/(JA) 



FIG. 6: The correlation function (5JS'J +m _ 1 )T divided by its zero-temperature value (5 f JS'J +m _ 1 )T=o 
for m = 2,3 in the massive regime A > 1. The dotted horizontal line gives the T = 0- value 1. The 
temperature is normalized with respect to J A. 
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